In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
from pathlib import Path
import os
In [3]:
import pandas as pd
import numpy as np
import sys
sys.path.append("../../")
DATA_PATH = Path('../../../data/')
MODEL_PATH = Path('../../../models/keras2lwtnn')
In [4]:
width = 6.693
height = 9.466
import matplotlib as mpl
mpl.use('pdf')
%matplotlib inline
import matplotlib.pyplot as plt
In [5]:
plt.style.use('seaborn')
In [6]:
plt.rc('font', family='serif', serif='Computer Modern Roman')
plt.rc('text', usetex=True)
# plt.rc('xtick', labelsize=8)
# plt.rc('ytick', labelsize=8)
# plt.rc('axes', labelsize=8)
# plt.rc('figure', titlesize=10)
In [7]:
from sklearn.model_selection import train_test_split
In [8]:
k = 1
data = pd.read_hdf(DATA_PATH / 'train_data.h5', 'train_set').sample(frac=k, random_state=137)
In [9]:
unused_features = [
'seed_nbIT', 'seed_nLayers', 'seed_mva_value', 'seed_nLHCbIDs',
'is_downstream_reconstructible_not_electron', 'is_true_seed',
'has_MCParticle_not_electron', 'has_MCParticle'
]
label_names = [
'is_downstream_reconstructible'
]
In [10]:
data.head()
Out[10]:
Training model for comparison¶
In [11]:
train_set, test_set = train_test_split(data, test_size=0.3, random_state=42)
In [12]:
x_train = train_set.drop(label_names + unused_features, axis=1)
y_train = train_set[label_names].copy().astype(np.int32)
x_test = test_set.drop(label_names + unused_features, axis=1)
y_test = test_set[label_names].copy().astype(np.int32)
In [13]:
def latex_rename(df):
d = {column:column.replace("_", "\_") for column in df.columns}
return df.rename(columns=d)
In [14]:
latex_rename(x_train).hist(bins=25, figsize=(width, height), layout=(5,2))
# plt.savefig("features_hist.pdf")
Out[14]:
In [15]:
temp_df = train_set.drop(label_names + ['seed_mva_value','is_downstream_reconstructible_not_electron', 'is_true_seed',
'has_MCParticle_not_electron', 'has_MCParticle'] , axis=1)
temp_df['log10_seed_p'] = np.log10(temp_df['seed_p'])
temp_df['log10_seed_pt'] = np.log10(temp_df['seed_pt'])
del temp_df['seed_pt']
del temp_df['seed_p']
latex_rename(temp_df).hist(bins=100, figsize=(width, height),layout=(5,2))
plt.tight_layout()
# plt.savefig("features_hist.pdf")
Random Forest¶
In [16]:
from sklearn.ensemble import RandomForestClassifier
forest = RandomForestClassifier(random_state=42, n_jobs=3, max_depth=20, class_weight='balanced')
forest.fit(x_train, y_train.values.ravel())
Out[16]:
In [17]:
from sklearn.metrics import roc_auc_score
print("train roc:{:.5}".format(roc_auc_score(y_train, forest.predict_proba(x_train)[:,1])))
print("train acc: {:.5}".format(forest.score(x_train, y_train)))
print("test roc:{:.5}".format(roc_auc_score(y_test, forest.predict_proba(x_test)[:,1])))
print("test acc: {:.5}".format(forest.score(x_test, y_test)))
DNN¶
Load model¶
In [18]:
architecture_file = (MODEL_PATH / 'lwtnn_v1_13' / 'architecture.json').as_posix()
weights_file = (MODEL_PATH / 'lwtnn_v1_13' / 'weights.h5').as_posix()
In [19]:
with open(architecture_file) as json_file:
json_string = json_file.readlines()[0]
In [20]:
from keras.models import model_from_json
In [21]:
model = model_from_json(json_string)
In [22]:
model.summary()
In [23]:
model.load_weights(weights_file)
Load new data¶
In [24]:
new_data = pd.read_hdf(DATA_PATH / 'test_data.h5', 'test_set')
new_data = new_data.drop(unused_features, axis=1)
In [25]:
x = new_data.drop(label_names, axis=1)
y = new_data[label_names].copy().astype(np.int32)
In [26]:
latex_rename(x).hist(bins=50, figsize=(width, height), layout=(5,2))
Out[26]:
In [27]:
weights = np.ones_like(y.values)/float(len(y.values))
plt.figure(figsize=(width/1.75,width/2.5))
plt.hist(y.values, weights=weights)
plt.title(latex_rename(y).columns[0])
# plt.savefig("class_balance.pdf")
Out[27]:
In [28]:
def nn_data_preprocessing(x):
x = x.copy()
x.loc[:, 'seed_angle'] = (np.arctan(x['seed_y'].values/x['seed_x'].values) + 0.00059458703227136336 ) * 1.3499114679506281
x.loc[:, 'seed_pr'] = (np.arctanh(x['seed_pt'].values/x['seed_p'].values) -0.2373458146766905) * 5.4252485447962835
x.loc[:, 'seed_r'] = (np.sqrt(x['seed_y'].values**2 + x['seed_x'].values**2) -712.28415274688871 ) * 0.0019285089504441804
x.loc[:, 'seed_chi2PerDoF'] = (np.sqrt(x['seed_chi2PerDoF'].values) -1.2984056703378601 ) * 2.3229264617966523
x.loc[:, 'seed_p'] = (np.log10(x['seed_p'].values) -3.8730802388444614) * 2.2146322148563526
x.loc[:, 'seed_pt'] = (np.log10(x['seed_pt'].values) -3.0783963469455355 ) * 4.6109558947469971
x.loc[:, 'seed_tx'] = (x['seed_tx'].values + 0.0024932173949750538) * 3.2600086889926385
x.loc[:, 'seed_ty'] = (x['seed_ty'].values + 0.00060514053765673619) * 14.811235511812143
x.loc[:, 'seed_x'] = (x['seed_x'].values -0.19295878936748806) * 0.0014091856183636067
x.loc[:, 'seed_y'] = (x['seed_y'].values + 4.7701536104670419) * 0.0019151935939470761
x = x[['seed_chi2PerDoF','seed_p','seed_pt','seed_x','seed_y','seed_tx','seed_ty','seed_angle','seed_pr','seed_r']].copy()
return x
x_nn = nn_data_preprocessing(x)
In [29]:
temp_df = x_nn.copy()
temp_df['log10_seed_p'] = temp_df['seed_p']
temp_df['log10_seed_pt'] = temp_df['seed_pt']
del temp_df['seed_pt']
del temp_df['seed_p']
latex_rename(temp_df).hist(bins=100, figsize=(width,height), layout=(5,2))
plt.tight_layout()
# plt.savefig("scaled_features.pdf")
In [30]:
x_nn.mean()
Out[30]:
In [31]:
y.hist()
Out[31]:
In [32]:
model.predict(x_nn.values[:5])
Out[32]:
Compare¶
In [33]:
from utils import plot_roc_curve, plot_confusion_matrix
from sklearn.metrics import roc_curve, confusion_matrix, accuracy_score
In [34]:
print("accuracy - Random Forest")
accuracy_score(y.values, forest.predict(x.values))
Out[34]:
In [35]:
print("accuracy - DNN")
accuracy_score(y.values, np.array(model.predict(x_nn.values) > 0.5, dtype=np.int32))
Out[35]:
In [36]:
predictions_df = pd.DataFrame(model.predict(x_nn.values), columns=y.columns, index=y.index)
In [37]:
plot_confusion_matrix(
confusion_matrix(y, forest.predict(x)),
classes=["ghost","true"],
title='Normalized confusion matrix - RandomForest',
normalize=True
)
In [38]:
plot_confusion_matrix(
confusion_matrix(y, np.array(predictions_df.values > 0.5, dtype=np.int32)),
classes=["ghost","true"],
title='Normalized confusion matrix - DNN',
normalize=True
)
In [39]:
fpr_forest, tpr_forest, thresholds_forest = roc_curve(y, forest.predict_proba(x)[:,1])
fpr_model, tpr_model, thresholds_model = roc_curve(y, predictions_df.values)
plt.figure(figsize=(8,8))
plot_roc_curve(fpr_forest, tpr_forest, "Random Forest")
plot_roc_curve(fpr_model, tpr_model, "DNN")
plt.legend(loc="lower right")
Out[39]:
In [40]:
# with open('roc.csv', 'w') as f:
# print('fpr,tpr', file=f)
# for i, (x,y1) in enumerate(zip(fpr_model, tpr_model)):
# if i%100 == 0:
# print("{},{}".format(x,y1), file=f)
In [47]:
import plotly.offline as py
py.init_notebook_mode()
import plotly.graph_objs as go
def plot_true_positives_and_negatives(
y_true,
probabilities,
normalize=False,
step=0.1,
title='True positives and true negatives vs threshold', ):
thresholds = np.arange(0.0, 1.0, step)
true_positives_rate = np.empty(thresholds.shape)
true_negatives_rate = np.empty(thresholds.shape)
for i, threshold in enumerate(thresholds):
classified_examples = np.array(
probabilities > threshold, dtype=int)
cm = confusion_matrix(y_true, classified_examples)
if normalize:
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
true_positives_rate[i] = cm[1, 1]
true_negatives_rate[i] = cm[0, 0]
plots = [
go.Scatter(x=thresholds, y=true_positives_rate, name='true positives'),
go.Scatter(x=thresholds, y=true_negatives_rate, name='true negatives'),
]
layout = go.Layout(
title=title,
xaxis=dict(title='threshold'),
)
fig = go.Figure(data=plots, layout=layout)
py.iplot(fig)
return true_positives_rate, true_negatives_rate, thresholds
In [54]:
plot_true_positives_and_negatives(
y, forest.predict_proba(x)[:,1],
title='Thresholds - RandomForest',
step=1e-2,
normalize=True
)
Out[54]:
In [55]:
true_positives_rate, true_negatives_rate, thresholds = plot_true_positives_and_negatives(
y, predictions_df.values,
title='Thresholds - DNN',
step=1e-2,
normalize=True
)
In [50]:
# with open('tpr_fpr.csv', 'w') as f:
# print('threshold,tpr,fpr', file=f)
# for x,y1,y2 in zip(thresholds, true_positives_rate, true_negatives_rate):
# print("{},{},{}".format(x,y1,y2), file=f)
In [51]:
print("ROC score Random Forest")
roc_auc_score(y.values, forest.predict_proba(x)[:,1])
Out[51]:
In [52]:
print("ROC score DNN")
roc_auc_score(y, predictions_df.values)
Out[52]:
In [53]:
cm = confusion_matrix(y, predictions_df.values > 0.5)
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
cm
In [ ]:
cm = confusion_matrix(y, predictions_df.values > 0.245)
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
cm
In [ ]:
from sklearn.metrics import log_loss
In [ ]:
log_loss(y, predictions_df.values)
In [ ]:
log_loss(y, forest.predict_proba(x.values)[:,1])
In [ ]:
dnn_pred = pd.DataFrame(model.predict(nn_data_preprocessing(x_train).values), columns=y_train.columns, index=y_train.index)
In [ ]:
dnn_pred['decision'] = dnn_pred['is_downstream_reconstructible'] > 0.5
In [ ]:
dnn_pred['correct'] = dnn_pred['decision'] == y_train['is_downstream_reconstructible']
In [ ]:
dnn_pred.head()
In [ ]:
def inv_expit(x):
return -np.log(1.0/x -1.0)
In [ ]:
inv_expit(dnn_pred[dnn_pred['correct'] == False]['is_downstream_reconstructible']).hist(bins=100, figsize=(width*1.5,width), alpha=.75, normed=True)
In [ ]:
np.max(inv_expit(dnn_pred[dnn_pred['correct'] == False]['is_downstream_reconstructible']))
In [ ]:
np.max(inv_expit(dnn_pred[dnn_pred['correct'] == True]['is_downstream_reconstructible']))
In [ ]:
np.min(inv_expit(dnn_pred[dnn_pred['correct'] == True]['is_downstream_reconstructible']))
In [ ]:
inv_expit(dnn_pred[dnn_pred['correct'] == True]['is_downstream_reconstructible']).hist(bins=100, alpha=0.75, normed=True)
In [ ]:
inv_expit(dnn_pred[dnn_pred['correct'] == False]['is_downstream_reconstructible']).hist(bins=100, figsize=(width*1.5,width), alpha=.75, normed=True)
inv_expit(dnn_pred[dnn_pred['correct'] == True]['is_downstream_reconstructible']).hist(bins=100, alpha=0.75, normed=True)
plt.title("Wrong decisions")
plt.rc('text', usetex=True)
plt.rc('xtick', labelsize=10)
plt.rc('ytick', labelsize=10)
plt.rc('axes', labelsize=10)
plt.rc('figure', titlesize=12)
plt.xlabel("$\hat{y}$")
plt.tight_layout()
# plt.savefig("wrong.pdf")
In [ ]:
dnn_pred[dnn_pred['correct'] == True]['is_downstream_reconstructible'].hist(bins=100, figsize=(width/2,width/2.5))
plt.title("Correct decisions")
plt.xlabel("$\hat{y}$")
plt.tight_layout()
# plt.savefig("correct.pdf")